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Abstract 



We present a method which, starting directly from QCD, permits a 
systematic gauge-invariant expansion to be made for all hard processes 
involving quarkonia in powers of the quark relative velocity, a small nat- 
ural parameter for heavy quark systems. Our treatment automatically 
introduces soft gluons in the expansion. Corrections arising from the in- 
corporation of gauge symmetry turn out to be important for decay and 
fragmentation processes involving QQ systems. The contribution of soft 
gluons is shown to be of higher order in v and so is neglected for calculations 
done upto and including 0{v^). 



INTRODUCTION 

The principles of quantum chromodynamics (QCD) were applied almost 20 years 
ago to the bound states of heavy quarks, such as cc and bb. These are pos- 
sibly the simplest strongly bound systems that exist. The large mass of the 
heavy quark sets a mass scale large enough so that perturbative QCD, together 
with a non-relativistic potential model description of the bound state, provides a 
good starting point to describe the decay and formation of quarkonia. However, 
quantitative predictions of the simple quarkonium model, even supplemented by 
radiative corrections, sometimes fail badly ^. Over the years, hundreds of papers 
^For a recent review of quarkonium phenomenology see, for example, Schuler^^). 
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have been written to rectify some of the failures. Nevertheless one still does not 
have a complete solution to this important problem of non-perturbative QCD. 

Our investigation into this venerable subject was prompted by the obser- 
vation that the fundamental symmetry to which QCD owes its origin, gauge 
symmetry, is manifestly violated by the naive quarkonium model. This is not 
hard to see: under a local gauge transformation q{x,t) — > U{x,t)q{x,t) the 
state normally used to describe quarkonia. 



does not remain invariant. In the above equation F is a space-time independent 
matrix in spin, color, and flavor indices and f{x) is the relative wavefunction. A 
gauge invariant state can be constructed, however, by inserting a path-ordered 
gauge link operator between quark operators. This amounts to including arbitrary 
numbers of soft gluons for transporting color between quarks. Building on this 
idea, in a previous publication*^^) we had proposed a manifestly gauge invariant 
effective field theory describing the interaction between heavy quarks, gluons, and 
quarkonium. The corrections accruing from the incorporation of gauge symmetry 
turned out to be substantial for decay and fragmentation processes, as well as 
radiative transitions, and this indicated the importance of a correct treatment. 
However, the relation of the effective theory to QCD was not transparent and it 
was not clear how the theory could be systematically extended to higher orders. 

In this paper we have developed a method which starts directly from QCD 




(1) 
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and which does allow for a systematic treatment of all high momentum transfer 
processes involving quarkonia such as inclusive decays, production, and fragmen- 
tation. The natural expansion parameter is the quark relative velocity which, for 
a heavy quark system, is small. All of the non-perturbative physics turns out to 
reside in a small number of matrix elements of gauge-invariant operators which 
are identified from symmetry considerations. The method used in this paper 
was inspired by the Feynman diagram treatment of deeply inelastic scattering as 
originally developed by Ellis, Furmanski, and Petronzio*^^-* , and recently further 
expanded upon by Jaffe and Ji*^^). 

While this work was nearing completion, we received a preprint authored by 
Bodwin, Braaten, and Lepage*^^) which presents a comprehensive QCD analysis of 
hard processes involving quarkonia. Their analysis is based upon a non-relativistic 
formulation of QCD. We share similar conclusions although these two approaches 
are quite different; with suitable identification of parameters the results are iden- 
tical. Perhaps an advantage of our method is its relative simplicity and its closer 
relation to the more familiar relativistic QCD. On the other hand, the work of 
Bodwin et al^^^ has wider scope because it is ultimately aimed at also generating 
the static properties of quarkonia through lattice calculations. We see the two 
approaches as complementary to each other. 
FORMALISM 

Our goal is to arrive at a systematically improvable, gauge-invariant, de- 
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scription of all hard processes involving a QQ system. By way of introduction, 
consider the decay of a positive C parity state into 2 photons (Fig. la) and the 
simplest Feynman graphs (Fig. lb) which contribute to it. Decay widths of spe- 
cific hadrons were computed from the leading approximation to these graphs long 
ago, and references may be found in Schuler*^^). For our purposes, it is useful to 
write the zero-gluon amplitude in Fig. lb as, 

= I ^JrM{k)h^''{k). (2) 

M{k) is the usual, but obviously non gauge- invariant, Bethe-Salpeter amplitude, 

M{k) = J d^xe^'=-^(0|r[V;(-x/2)V'(x/2)]|P). (3) 

The tensor h^^{k) is the amplitude for two quarks, not necessarily on their 
mass-shells, to annihilate into 2 photons. To leading order in as this is, 

h^'•'{k) = -ie^ [YSpik + K)'y^' + -f^Spik - K)Y] ■ (4) 

In Eqs 2-4, x'^ is the relative distance between quarks, k'^ — |(pi — P2Y 
quark relative momentum, K^^ = ^{q2 — Qi)^ is the photon relative momentum 
and Sf is the free fermion propagator. We shall refer to M{k) and h^''' as "soft" 
and "hard" parts respectively in the following. 

Now consider the fact that a large momentum of 0{m), where m is the 
heavy quark mass, flows through the single propagator in Fig. lb but that, on the 
other hand, the soft part has typical quark momenta much less than m. This 
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suggests that we expand the hard part in powers of k°', 



h'''{k) = Y.^k-\..k--UZ...^„. (5) 



where, 

d _d 



UZ..a. = ^■■■■^hnk=o. (6) 



Inserting Eq.5 into Eq.2 and integrating by parts gives for the amphtude, 

U''{k) = Tr^M''-''"UZ..ar., (7) 

n 

where, 

1 _ <-»ai <->an 
^a,...ar. = — 0|^ id .... id 1p\P). (8) 

n\ 

The matrix elements in Eq.8 have derivatives d = ^{d —d ) evaluated at zero 
relative quark separation. As we shall see later, in a non-relativistic model the 
n = matrix element is proportional to the wavefunction at the origin, and so 
on. However, we do not need to appeal to any particular model at this point. 

Next consider the single gluon diagram in Fig. 2. The corresponding ampli- 
tude is, 

^ / (0^ ^')^'''''(^, ^0- (9) 

The "soft" part Mp{k, k') is a generahzed B-S amplitude. 



Mp{k,k')= J d^xd^z e"'-''e"'-'{0\T[iP{-x/2)Ap{z)ilj{x/2)]\P), (10) 

and Ap = ^X"'A'^ is the gluon field matrix. The "hard" part H^'^'^ is the annihi- 
lation amplitude for QQg — > 77. To leading order this is, 

6 



H^'''f'{k, k') = -ie^g WSpik + h' + K)YSF{k - \k' + K)^^' 

Hi^^v). (11) 

Excluded are the diagrams such as in Fig. 2b. These are properly included in 
Fig. lb since the lines emerging from the blob correspond to interacting fields. 
Expanding the hard part, 

H'^^'^ik, k') = (12) 

where, 



(9 3 3 3 



An integration by parts on k,k' yields an alternate form for the amplitude , 
Tr{k) ^TtY: M----'^-^'K'Za„A...A' (14) 

nl 

where, 

^ai...a„,/3i...A = -|T7(0|^ i 3 ....i 3 ipi3'^' 3^'Ap\P). (15) 

The derivatives i 3 act only upon the quark operators. 

Finally consider the two-gluon contribution to the amplitude shown in Fig.3. 
The amplitude corresponding to Fig.3a is, 

= / (16) 

where, 

{0\^{-x/2)A,,{x')A,.{x")^|;{x/2)\P). (17) 



The product of color matrix fields is appropriately symmetrised in the above 
because of Bose symmetry. The hard part H^^^^'f" is, at lowest order, 

H^^^'^XK k', k'')^-ie^g^[rSF{k+^''+Ky''SF{k+k'-^"+Ky'SF{k- ^k"+K)^'' 

+crossed]. (18) 

This may be expanded as before about k — k' — k" — 0. 

We now make the observation that the simple Ward identity, 

^ Sf{p) = -Spiph^Spip), (19) 



Qpa 

leads to a number of useful relations. In particular, 

V^"" = -gU^'^f, (20) 

allows us to combine the n = 1 term in Eq.7 and the n = I = term in Eq.l4 
into a gauge invariant sum, 

M'^U^" + M'^V^'^ = (0|^ iDi^lPW^". (21) 

Similarly, the leading order term in the hard two-gluon amplitude is just the 
second order term in the zero-gluon amplitude, 

^m^pV'(o) (22) 

Collecting together appropriate terms leads to another gauge-invariant matrix 
element, 

M^-'UZ' + M^Vr' + Mp,p.Hr''"{0) = l^iDp' iDp'^W'"''''" . (23) 



Here Mptpt> is just the leading order term from Eq.l7. It is the matrix element 
with all fields at the same space-time point, 

M^V" = (OIV^A^'WI^)- (24) 
Next, look at the n — Q I — 1 term in Eq.l4: 

{Q\^l)il]idpAp\P)d'f^H''^p = -{Q\il)il]{dpAp - dpAp)\P)d'^m^P (25) 

The last step made use of d'^H^^^P — —d'^H^"^ , with all derivatives evaluated at 
k = k' = 0. The quantity inside the brackets is the abelian field strength tensor; 
the non-abelian part, ig [A13, Ap], can be shown after some effort to arise from the 
amplitude in Fig. 3b. 

We now collect results together and summarize. The contribution from 
the 0, 1, 2 gluon diagrams in Figs. 1-3 have been expanded in powers of relative 
momentum, added together, and terms suitably rearranged. The total amplitude 
is, neglecting higher powers of momentum. 



(To + Ti + Ts)'^" = Tr 



{0\ipip\P)h'''' + {0\ipi Da ip\P)d'^h'"' 
+{0\ipiDa iD^ip\P)-d''d^h^"' 



(26) 



The hard amplitude hP^{k), H^^^^^k, k'), and their derivatives are all evaluated at 
k — k' — 0. We see that each term in the above is a product of a gauge invariant 
matrix element characteristic of the decaying hadron and a simple, calculable. 



hard part. In the following, hadrons of specific J^^ will be considered and the 
relative order of importance of the terms in Eq.26 will be explicated. Radiative 
corrections, which are not included in the lowest order amplitudes h and H, will 
be considered separately. 



MATRIX ELEMENTS 

In the previous section the amplitude for a C—-\- quarkonium state to de- 
cay into 2 photons was expressed in terms of matrix elements of leading gauge- 
invariant operators. Further progress requires we specify the angular momentum 
and parity: we take J=0, V—— {rjc and r)i, mesons) for now, leaving other mesons 
for later analysis. From Lorentz invariance, and invariance under charge conju- 
gation and parity, the only non-zero matrix elements are, 

(0|^^|0-+) = aiM275 + asMfys (27) 

(0|VSiB?^|0-+) = ifcMV^'^TsP^ (28) 

(0|V'F''''V|0"+) = cMh'"'''>^-faP/3 (29) 

(0|V^iBzB^|0-+) = M^\diM^g'"'+d2P^P'']'y5 + 



pa p/x pv 

e.g^'P'^+e ^ +e3(^"^P'-+g"'-P^)^a75 (30) 



For brevity, color has not been explicitly indicated in Eqs 27-30. It is clear that 
in Eq.27 the two quarks must be in a color singlet and so, regarded as matrix 
in color space, only the unit operator appears on the right hand side. However, 
when one gluon appears, as in Eqs 28-29, the quarks may be in either singlet or 
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octet states and the corresponding constants b^^\..b''^^ then appear on the RHS. 
In Eq.30, the two gluons can combine into a color singlet or octet, and those in 
turn can combine with the quark singlet and octet respectively to give an overall 
singlet. Clearly this leads to a large number of constants which must be known 
in order to describe corrections to 0~+ decay and if our approach is to have any 
practical utility, this number must be curtailed according to some well defined 
principle. 

To progress beyond this point, it will be necessary to specialize our hitherto 
general discussion and select a particular gauge. The Coulomb gauge is natural 
for this problem, as shown by vast experience with positronium states. We shall 
not repeat here the arguments of Lepage et al*^^) who, using the QCD equations 
of motion in the Coulomb gauge, make the following estimates, 

— * ~* 9 

d ~ mv , gA ~ mv 
gE^m\^ , gBr^m^v^. (31) 

Here v is the relative velocity of quarks - the small parameter in the theory. The 
estimates (31) allow us to see that explicit gluons will not enter in the leading 
order corrections to the naive quarkonium model. Therefore, working to 0{v^), 
one may effectively replace the covariant derivatives in Eqs 27-30 with ordinary 

— * — * 

derivatives, and ignore E and B. 

We next observe that tracing Eq.30 with 757^ or 757,^, and using the equa- 

11 



tion of motion ilfkl) — mip, yields the constraint, 



€1 + 62 + 5es = 0. 



(32) 



Working in the rest frame of the meson P^^ — (M, 0) and putting /j, — v — i 



yields 63 ~ 0{v^). Hence ei = —62 + O(v^). Setting fj, — 1/ — yields d2 — 
—di + 0{v^). This leaves us with having to deal with 01,02,6, di and Ci - five 
independent parameters at the O(v^) level. 

Further progress demands that we specialize a step further and specify a 
model framework for the O""*" quarkonium state. We shall assume, in common 
with many other authors, that the Bethe-Salpeter equation with an instantaneous 
kernel does provide an adequate description. This has been conveniently reviewed 
by Keung and Muzinich^^^ and we adopt their notation. The momentum space 
B-S amplitude x{p) satisfies the homogeneous equation. 



which, after making the instantaneous approximation K{P,p,p') — V{p,p') and 
reduction to the non-relativistic limit yields^*^^ , 




(33) 



x{p) 



M^^jM -2E){E + m- ^.7)75(1 --fo){E + m- p.7)0(|p1) 
4:E{E + m){pO + f-E + ie){p^-f + E- ie) 



(34) 



The scalar wavefunction 4>{\p\) is normalized to unity. 




(35) 



and. 




(36) 
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Fourier transforming x{p) to position space yields (O|'0(— x/2)'0(x/2)|P) from 
which, by tracing with appropriate gamma matrices, the coefficients oi, 02,^,^1 
and ei can be extracted. So finally, to O(f^), one has a rather simple result, 

(OI^^IO--^) = iMVV(0)(l + |)75 + M-V^^M^fy5, (37) 

(0|^^Z^^^|0-+) = iMV2Z!M,^M.^^p^, (38) 

(0|i/;F'^>|0-+) = 0, (39) 

mD^iD'^^lO-) = ^^^/^^(^'^^-^) (1 + 1)7. (40) 

Eqs. 37-40 express hadronic matrix elements, upto 0{v'^), in terms of two basic 
parameters: 0(0) and V^0(O). These may be obtained for any given phenomeno- 
logical potential from a non-local Schrodinger type of equation^^^ . 

DECAY RATES 

The decay rate for ^ — > 27 may be directly computed from Eqs. 26 
and 37-40. The calculation is facilitated by the observation that, from invariance 
under time reversal, the crossed diagrams in Figs. 1-3 exactly double the uncrossed 
ones. The result is, 

+ ^' + - M3/^(i^+^^) (*(°) + 

The factor of \/3 comes from the sum over colors. The quark mass m differs from 
M/2 because of the strong binding, 

= 2m - M. (42) 
13 



sb/M is of 0{v'^) from virial theorem, and thus of the same order of magnitude 
as V^/M^. From Eq.41 it is simple to get the decay rate (excluding radiative 

corrections) , 

V^ + Vb + Tc + Tr. (43) 

In Eq.43, Vq is the conventional result, 

r„ = ^fl^(O), (44) 

where eg is the quark charge and i?(0) = 0(O)-\/47r. is the correction coming 
from m ^ M/2, 

= -2^Fo, (45) 

and Tc is the term coming from differentiating the quark propagator once, and 
then twice, 

16 V^i?(0) 

^^-3M^"R(0r ^ ^ 

Lowest order radiative corrections to — > 27 were calculated by Barbieri et 

al^*'' a long time ago. These are 0{v'^) too. 

For decay into 2 gluons, the only difference in Eqs. 44-46 is from the color factor, 
but the 3 gluon vertex changes the form of the radiative correction, 

ro-+-.2, = (Fo + Fb + Fc + F^) , (48) 

where, 

F', = log £ + 1^ - - 11 log 2 + n.il log 2 - ^)) ^^Fo- (49) 
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The radiative corrections to the decay into gluons involves both the renormal- 
ization scale /i and the renormalization scheme; for a discussion of this point see 

Kwong ct al . 

1— DECAY 

The formalism developed for two photon decay can be used quite trivially 
to calculate the important decay 1 — 7* — The "hard part" is 
the single, momentum independent vertex, h'^ — —ie^^. There are therefore no 
corrections from expanding the hard part, and the amplitude analogous to Eq.26 
is simply, 

T^" = Tr{Q\^l:il:\P,e)h^, (50) 

where e'^ is the vector meson polarization vector. Going to the Coulomb gauge, 
and reducing the B-S equation, yields the amphtude analogous to Eq.34 with the 
simple replacement 75 — > Using C and V invariance of the matrix element in 
Eq.l, we find that to 0{v'^), 

(Ol^^lP.e) = iMV^d + Zl)^(0)(l + L) ^_ iMV^^d - ^) ^ (51) 
This yields for the decay to leptons, 

Tl — -^l+l- — ^VW + Trad + Tcor- (52) 
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FyvK is the usual Van Royen-Weisskopf^-^^^ formula^, 

^vw-^R^O), (53) 

Trad is the radiative correction calculated some time ago by Celmaster^^^^ , 

16 

Trad = —-;;—Ois^vw, (54) 

OTT 

and Fcor is the correction term which comes from Eqs.50 and 51, 

Although we have used the same symbol R{r) for the radial wavefunction of 
the 1 and 0~+ states, these wavefunctions are in principle different. We shall 
return to this point later. 



COMPARISON 

In two important previous works, O(v^) corrections to and 1 quarko- 
nium decays have been evaluated. The first approach by Keung and Muzinich^^^ 
starts from the B-S equation with an instantaneous kernel . Subsequently a non- 
relativistic reduction is made, followed by an expansion of the lowest order am- 
plitude about the mass-shell value of the relative momentum p ^ = (M/2)^— = 
—m€B- The relevant results of their work are in Table 1. Their treatment does not 
satisfactorily resolve the issue of QCD gauge invariance of decay rates, although 
they do raise this question. 

^The lepton mass correction is simply included by multiplying Eq.53 by — -^^(1 + 1^) 
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The second approach is that of Bodwin et al^^^ which builds systemati- 
cally upon the rigorous formulation of non-relativistic QCD by Lepage and co- 
workers*^^-*. These authors introduce an ultra-violet cutoff A of 0{m) and then 
construct a NRQCD Lagrangian by successively adding new local interactions 
with 2-component fermion spinors. To account for annihilation into photons, 
higher dimensional terms involving 4 fermion operators are introduced into the 
Lagrangian, and their coefficients are computed in a power series in by putting 
the annihilating quarks on mass-shell. The annihilation process, which cannot be 
got directly from NRQCD, is taken into account via the optical theorem which 
relates annihilation rates to the imaginary parts of QQ — > QQ scattering am- 
phtudes. Bodwin et al^^^ express their results (see Table 1) in terms of non- 
relativistic wavef unctions, their derivatives and the quark mass m. They do not 
use the meson mass M. However, to enable a comparison, we have expressed 
their results in terms of M using €3 — 2m — M after expanding to first order in 
eb/M. Note also that this definition of is opposite in sign to that of Keung 
and Muzinich^^) . 

The third approach is that of this paper. For completeness we summarize 
this too: the decay amplitude is given by the sum of all distinct Feynman diagrams 
leading from the initial quarkonium state to the final state. Each diagram is put 
into the form of a (multiple) loop integral with a kernel which is a product of a 
hard part and a soft part. The hard part is treated with perturbative QCD, and 
the soft part is analyzed into its different components with the use of Lorentz, 
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C, and V symmetries. Use of the QCD equations of motion enables separation 
of these components according to their importance in powers of v. As the last 
step, a specific commitment to dynamics is made and the B-S equation is used 
to express the components in the form of wavef unctions. 

The first comment regarding the results summarized in Table 1 is that all 
six entries collapse into a single one, 1 + ^j^, upon making the identification^ 
es/M — 2V'^R/M'^R. It is interesting to note that this condition is precisely that 
which follows for a potential V{r) which has V^(0) = 0. The Schrodinger equation 

for this potential at r = is , 

- |:V^i? = -sbR. (56) 

However, it is well known that at small r the potential is Coulomb-like, 

V{r) ~ 1/r. In this case the entries in Table 1 are not identical for arbitrary 

choices of or equivalently, the quark mass m. Furthermore, V^R is apparently 

singular at the origin V^-R(r) ~ MR{r)/r. As is clear from the uncertainty 

principle, the local kinetic energy becomes very large at short distances and the 

expansion in powers of v breaks down. This difficulty may be circumvented by 

imagining that annihilation takes place in a diffused region of size 0(l/m), i.e., 
^The relation between the binding energy and is explained briefly as a renormalization 
condition in Labelle et al^^^^ in the NRQED approach (see their Eqs.ll and 12). However, in 
our case there is no principle which a priori constrains es to bear a fixed relation to and 
therefore both will be considered adjustable parameters. 
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that R and V^-R are quantities renormalized at this scale. In any case, they are 
simply parameters which serve instead of the parameters in Eqs. 27-30. 

In order to estimate the correction factors for charmonium, we used the 
following values of the independent parameters, 

as = 0.19 
m = 1.43 GeV 
^ = -0.7 GeV. 

With this particular choice of parameters and using Eqs. [43], [48] and [52], the 
decay rates are calculated to be, 

T{J/ilj — >e+e-) = 5.61 KeV 
r{r]c — > hadrons) = 9.99 MeV 

r(77c — > 27) = 6.48 KeV. (57) 

These values agree reasonably well with the experimentally measured decay widths 
which are*^^^'', 

r{J/jjj — > e+e") = 5.36 ± 0.28 KeV 
r{r]c — > hadrons) = 10.3 ± 3.6 MeV 

r{ric — ^27) = 8.1 ±2.0 KeV. (58) 

In evaluating expressions (57), the radiative corrections are calculated at the 
renormalization point /i — m^^\ The wavefunctions of J/%1) and r]c at the origin 
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differ from each other to 0{v'^). This difference is neglected in taking the ratio 
Their values are, 

= 0.978 GeV^ 
\Rr,f ^ 0.936 GeV^. 

One remark concerns the value of as used above, which differs from the 
value deduced from deep inelastic scattering, as{mc) ~ 0.3. The reason is the 
following: the value of the parameter depends upon the value of chosen 
and, for smaller values of ag, this is negative. The corresponding values of the 
wavefunctions of J/ip and rjc at the origin should differ from each other by 
by the assumptions used in this paper. However, for larger values of q;^, 
becomes positive and correspondingly the difference between the wavefunctions 
becomes rather large. For example, for q;s=0.24, we have, ^^ = 2.8 GeV^, with 
I i?^/^ 1 2= 0.582 GeV^ and |i?^J2=0.194 GeV^. The difficulty in using large values 
of as has also been noted by Consoli and Field and suggests that 0{al) 
radiative corrections to charmonium decays may well be significant. 

In conclusion, we have investigated higher order corrections to the decay of 
O""*" and 1 heavy quarkonia and shown how these corrections can be system- 
atically incorporated in terms of various bound state matrix elements of gauge- 
invariant quark and gluon operators. Investigation of V—+, C—+ states - which 
correspond to P-waves in the n.r. limit - is in progress. We are also currently 
calculating, using the framework developed in this paper, the more complicated 
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case of the decay of negative C parity quarkonium into 3 gluons/photons. This 
will enable a more detailed comparison of theory vs experiment. 
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Figure Captions 

Figure:! a) Heavy quarkonium decay into two photons, b) The lowest order 
diagrams contributing to QQ — > 2j. 

Figure: 2 a) Single gluon diagrams, b) Interacting quark field diagram which is 
properly included in Fig. lb. 

Figure: 3 a) Two gluon diagrams, b) Three gluon vertex diagrams. 
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0-+ — > 2g 


I— —r j+r 


Keung et al^^^ 


^ ^ 3M 


^ 3M 


Bodwin et al^^^ 


1 2eB 1 16 V^i? 
^ M ^ 3M^ R 


1 2eB 1 16 V^i? 
M ^ 3M^ R 


This work 


1 2eB 1 16 V^i? 
^ M ^ 3M^ R 


^ ^ 3M^ 



Table 1: A comparison of the O(f^) correction factor, excluding radiative correc- 
tions, which multiply the zeroth order formulae for the electromagnetic decay of 
quarkonium states. M is the hadron mass, = 2mr-M is the binding energy, and 
both R and V^i? are evaluated atr — 0. Note that all six entries become identical 
upon making the identification eb/M — 2V^R/ M'^R. 
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